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ABSTRACT 

We test the ability of the numerical action method (NAM) to recover the indi¬ 
vidual orbit histories of mass tracers in an expanding universe in a region of radius 
26 /i“^Mpc, given the masses and redshift-space coordinates at the present epoch. 
The mass tracers are represented by dark matter haloes identified in a high resolution 
A-body simulation of the standard ACDM cosmology. Since previous tests of NAM at 
this scale have traced the underlying distribution of dark matter particles rather than 
extended haloes, our study offers an assessment of the accuracy of NAM in a scenario 
which more closely approximates the complex dynamics of actual galaxy haloes. We 
show that NAM can recover present-day halo distances with typical errors of less than 
3 per cent, compared to 5 per cent errors assuming Hubble flow distances. The total 
halo mass and the linear bias were both found to be constained at the 50 per cent level. 
The accuracy of individual orbit reconstructions was limited by the inability of NAM, 
in some instances, to correctly model the positions of haloes at early times solely on 
the basis of the redshifts, angular positions, and masses of the haloes at the present 
epoch. Improvements in the quality of NAM reconstructions may be possible using 
the present-day three-dimensional halo velocities and distances to further constrain 
the dynamics. This velocity data is expected to become available for nearby galaxies 
in the coming generations of observations by SIM and GAIA. 
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1 INTRODUCTION 


a nonlinear process, making numerical methods an essential 
tool for understanding the observed large-scale structure. * 


In the standard cosmological paradigm, the formation of 
large-scale structure is driven by the gravitational amplifica¬ 
tion of small initial density fluctuations (e.g. Peebles 1980). 
In addition to gravity, hydrodynamical processes can influ¬ 
ence the formation and evolution of galaxies, groups and 
clusters of galaxies. But since hydrodynamical effects play a 
minor role on scales larger than the size of galaxy clusters, 
gravitational instability theory alone can directly relate the 
present day large-scale structure to the initial density field 
and provide the framework within which the observations 
can be analyzed and interpreted. Gravitational instability is 


There are two complementary numerical approaches 
to studying cosmological structure. The first relies on N- 
body techniques designed to solve an initial value problem 
in which the evolution of a self-gravitating system of mas¬ 
sive particles is determined by forward numerical integration 
of the Newtonian differential equations. Because the exact 
initial conditions are unknown, comparisons between these 
simulations and observations are mainly concerned with gen¬ 
eral statistical properties. 


The second approach works in the opposite direction, 
deriving from the observed present-day distribution and pe- 


* Present address : Racah Institute of Physics, The Hebrew Uni¬ 
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culiar motions of galaxies, and independently of the nature 
of the dark matter, certain features of the dynamics at earlier 
times. The numerical action method (NAM) belongs to this 
second category of approaches. It arises from the observation 
that the present-day distribution of galaxies, combined with 
the reasonable assumption that their peculiar velocities van¬ 
ish at early times, presents a boundary value problem that 
naturally lends itself to an application of Hamilton’s prin¬ 
ciple in which stationary variations of the action are found 
subject to the boundary conditions. The result is a predic¬ 
tion of the full orbit histories of individual galaxies, either 
with real space boundary conditions (Peebles 1989, 1990, 
1994, 1995) or, after a coordinate transformation, in red- 
shift space (Peebles et al. 2001, Phelps 2002). 

The potential of NAM as a probe of galaxy dynamics 
and of cosmological parameters has been explored in a num¬ 
ber of studies following the introduction of the method in 
Peebles 1989. Possible applications include the full nonlinear 
analysis of orbit histories of nearby galaxies (Peebles 1990, 
1994, 1995; Sharpe et al. 2001), recovering the initial power 
spectrum of density fluctuations (Peebles 1996), predicting 
the values of cosmological parameters (Shaya et al. 1995), 
and estimating the proper motions of nearby galaxies (Pee¬ 
bles et al. 2000). Concerning the latter application, ground 
and space-based observations will soon make possible the 
measurement of the full three-dimensional velocities of many 
nearby galaxies and promise both a rigourous test of NAM 
predictions and, given the additional dynamical constraints 
on galaxy motions, the possiblity of using NAM as a probe 
of individual masses of nearby galaxies. 

Since a central result of NAM, the past orbit histories 
of galaxies, cannot be confirmed by direct observations, N- 
body simulations provide an important test of NAM and its 
key assumption that galaxies can be approximated as dis¬ 
crete, non-merging objects throughout their history. It is de¬ 
sirable then to test NAM in a scenario which approximates 
the complexity of the observational situation but where all 
of the relevant physical quantities are known. Previous tests 
of NAM using A^-body simulations have either been con¬ 
fined to a few dark matter haloes at the scale of the Lo¬ 
cal Group (Branchini & Carlberg 1994, Dunn & Laflamme 
1995), traced the paths of individual dark matter particles 
rather than extended haloes (Nusser & Branchini 2000), or 
used simulations which demonstrate in principal the ability 
of NAM to recover particle orbits to a high degree of ac¬ 
curacy but which do not reproduce the full complexity of 
extended mass distributions (Phelps 2002). 

In this paper we extend the tests of NAM to simula¬ 
tions at a scale approaching that of the local supercluster 
with a catalogue containing several hundred extended ob¬ 
jects modelled as particles. We begin with an overview of 
the relevant properties of the A^-body simulation and the 
halo catalogue we derived from it, and follow with details of 
the version of NAM used here, which includes a novel ap¬ 
proach to the assignment of halo masses. We will then test 
the sensitivity of NAM both as a probe of the total mass 
as well as of the linear bias, and examine in some detail a 
representative solution, focusing on the comparison between 
the NAM predictions and the actual halo orbits. 



r / h 'Mpc 


Figure 1. The average matter density f2(< r) in a sphere of 
radius r from the center of mass of the simulation, as a function 
of redshift. The arrows mark the nominal initial boundary (r = 
26 /s“^Mpc) of the high-resolution region. The average density 
within 26 /i“^Mpc is about 10 per cent larger than the cosmic 
mean flm = 0.3 at z = 0. 


2 THE SIMULATION 

We will test NAM in a A-body simulation of a ACDM cos¬ 
mology with matter density = 0.3, vacuum energy den¬ 
sity Ha = 0.7 (in units of the critical density pc), and a 
Hubble constant h = 0.7 (in units of 100 kms“^Mpc“^). 
The simulation has an initial fluctuation power spectrum of 
spectral index n = 1 and a present-day normalization am¬ 
plitude as = 0.9. It was run with the GADGET A-body 
code (Springel, Yoshida & White 2001). It is part of a suite 
of simulations which zoom in on a spherical region of ini¬ 
tial comoving radius 26 /i~^Mpc selected from a larger par¬ 
ent simulation (see Stoehr 2003 for details). This “zooming” 
technique has the advantage of preserving the large-scale 
power of the parent simulation while allowing the inner re¬ 
gion to be resolved at smaller scales. In particular, the mass 
of the 168436 “high-resolution” particles which sample the 
inner region is = 6.82 x 10^° Mq/Ii, in contrast to 

the 602272 “low-resolution” particles (each having a mass 
^ tUhr) which sample the outer region. 

Fig. 1 shows the average matter density H(< r) aver¬ 
aged over a sphere of radius r as a function of redshift. At 
the initial time step (z = 20), the mean density in a sphere 
of radius r = 26 /i^'^Mpc, H(< 26), is very close to Hm, but 
increases by about 10 per cent from initial to final time. The 
net mass inflow across the region boundary has the poten¬ 
tial to disrupt the reconstruction at the edge of the region 
and provides a further opportunity to test the sensitivity of 
NAM to the total mass. 
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M / h-'Mg 

Figure 2. The friends-of-friends differential mass function of dark 
matter haloes (in logarithmic bins Alogj^gm = 0.2) computed 
from the snapshot at 2 = 0 (histogram) compared to the theo¬ 
retical prediction of Sheth & Tormen (dashed curve). The arrow 
indicates M* ~ 10^^ M-q/K, or the typical mass of haloes collaps¬ 
ing at the present time. 

2.1 The halo catalogue 

The principal objective of this paper is to test the ability 
of NAM to reconstruct the orbits of simulated dark matter 
(DM) haloes which resemble the galaxy -|- halo distributions 
of the “real world”. To this end we need to specify the po¬ 
sitions, velocities and masses of a fixed number of haloes as 
a function of time, since our implementation of NAM con¬ 
serves particle number. However, haloes do not form at any 
single redshift, and mergers potentially play an important 
role in the dynamics. 

In order to create a catalogue of a fixed number of haloes 
from the inner region of the simulation, we first identify all 
the DM haloes above a given mass threshold Mmin in the 
final simulation time step 2 = 0, using a friends-of-friends 
(FOF) group-finding algorithm based on the method given 
in Davis et al. (1985). Only groups containing at least 20 
DM particles are classified as haloes, resulting in a mass 
threshold Mmin = 1-36 x 10^^ Mq//i. We used the stan¬ 
dard linking length of 20 per cent of the mean interparticle 
distance (Jenkins et al. 2001). Using haloes defined with 
a linking length of 30 per cent we found a significant im¬ 
provement in the accuracy of NAM reconstructions, which 
is not unexpected since this probes a more linear regime. In 
this paper however we confine our analysis to the standard 
linking length since it most closely reproduces the expected 
mass function and thus gives the best approximation to real 
galaxy haloes. 

The center of mass (CM) positions and velocities of all 
those DM particles which are associated with haloes at 2 = 0 
define the resulting catalogue of 576 objects which is the in¬ 


put to the numerical action code. The orbit histories of these 
objects are then defined as the CM positions and velocities 
of the DM particles which belong to the corresponding halo 
at 2 = 0. For our purposes the redshift outputs at 2=20, 4, 
3, 2, 1 and 0 were sufficient for comparison of halo orbits 
with the action reconstruction. 

To model the tidal field from the outer region of the 
simulation, we created an auxiliary particle catalogue from 
a random sample of 1000 low-resolution DM particles, multi¬ 
plying the masses by the appropriate factor of 277.9 to bring 
the mass density in the outer region to Dm = -3. To reduce 
the computational burden, the positions of these tidal field 
particles were fixed in the NAM reconstructions and their 
masses were evolved according to linear perturbation theory, 
as described in Shaya et al. (1995). Parallel sets of NAM so¬ 
lutions were created with and without this tidal field. 

The final step is to identify one of the haloes near the 
center of the simulation as the reference halo, and to build 
an input catalogue to NAM using the angular positions and 
redshifts relative to this halo. The analysis below was re¬ 
peated for two haloes within 8 Mpc/h of the CM of the 
high-resolution region, as a check on the sensitivity to the 
position and velocity of the reference halo, and since the re¬ 
sults were consistent we will confine our analysis to a single 
reference halo. 

In Fig. 2 we show as a histogram the FOF differential 
mass function of DM haloes in the high-resolution region at 
2 = 0. The arrow marks M* ~ Mq/Zi, the typical mass 
of haloes that have just formed by redshift 2 = 0. As can be 
seen, the halo masses span more than two orders of magni¬ 
tude, with the largest haloes having a mass M Mq//i. 

The dashed curve is the theoretical prediction of Sheth & 
Tormen (1999) for a ACDM cosmology matching that of 
the simulation. The good agreement shows that the high- 
resolution region of the simulation is a fair reproduction of 
the expected mass function. 

In the top panel of Fig. 3 we compare the auto¬ 
correlation of DM particles (solid) and haloes (dashed), 
and respectively. In the lower panel we plot the bias 
= Chh/^mm- The dotted horizontal line marks 6=1. 
Results are shown at redshift 2 = 0 (square) and 1 (tri¬ 
angle). Note that, on the scale r ^2 6“^Mpc, DM haloes 
are as biased as the matter, i.e. 6 ~ 1. This follows from 
the fact that at 2 = 0 most of the haloes have a mass 
M « lO’^^ Mq/6 (see Fig. 2), where Mi,{z) is the 

typical mass of haloes which collapse at redshift 2 . The bias 
factor inferred from the unsmoothed and character¬ 
izes the clustering of DM haloes. We will propose below that 
it is possible using NAM to directly measure this bias factor 
when the total mass is known. 


3 THE ACTION RECONSTRUCTION 

The implementation of NAM used in this paper is fully de¬ 
scribed in Peebles et al. (2001) and Phelps (2002). We refer 
the reader to appendix A for details. NAM solves for the or¬ 
bits of mass tracers given the boundary conditions that the 
object positions in redshift space are fixed at the present 
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Figure 3. The auto-correlation of haloes (dashed curve) and dark 
matter (solid curve) obtained from the unsmoothed halo and mat¬ 
ter distribution at ^ = 0 (square) and z = 2 (triangle). The bias 
is plotted in the lower panel. The horizontal dashed 
line shows 6 = 1. 


time and that at the initial time a^midx.i/dt —> 0. While 
some implementations of NAM expand the orbits in terms 
of suitable trial functions (Peebles 1989, 1990, 1994; Nusser 
& Branchini 2000), the present implementation relies on ma¬ 
trix inversion of a discrete form of the equations of motion 
(Peebles 1995, Phelps 2002), which allows the creation of 
more complex orbits. In this formulation the initial veloci¬ 
ties ad’Kijdt are not mathematically constrained to vanish. 
Our approach to this problem will be discussed in the fol¬ 
lowing section. 

The generation of solutions begins with a set of random¬ 
ized halo orbits with 10 time steps and relaxes to a target 
in the gradient of the action many orders of magnitude be¬ 
low that of the initial guess by iterated per-particle matrix 
inversion. Different choices of initial randomized orbits may 
yield different solutions, since uniqueness is not guaranteed, 
and so for each choice of cosmological parameters a suite of 8 
solutions was generated and the solution with the lowest 
(see below) was used. Orbits with negative radial distances 
are allowed mathematically but in such cases the orbit is 
re-randomized and the matrix inversion procedure repeated 
until all of the particles have positive distances. The result¬ 
ing orbit paths satisfy the standard equations of motion, as 
confirmed by a simple leapfrog integration forward in time 
from the first time step of the solution. Solutions with inter¬ 
secting orbits, which are also mathematically allowed, are ef¬ 
fectively suppressed by modelling each particle as a constant 
density sphere with radius proportional to the cube root of 
the mass, normalized so that a halo of mass 2.4 x 10 ^^Mq 
has a radius of 180 kpc. 

To more closely parallel the observational situation, 
where redshifts are known to greater accuracy than the dis¬ 


tances, the standard action was recast with a partial trans¬ 
formation of coordinates, as described in Phelps 2002, so 
that the redshift at a = 0 is fixed as the boundary condition 
instead of the radial distance. The distances then emerge as 
predictions, from which a standard measure is defined to 
measure the goodness of fit of the solutions: 


E / mod . cat\2 

(Mi - Mi ) 


( 1 ) 


where a-nd are the predicted and catalogue dis¬ 

tance moduli for each particle and Oi is the dispersion in 
the observed distance moduli. Since there are no observa¬ 
tional errors in the simulation, we added to the input dis¬ 
tance moduli a gaussian random error with a dispersion of 
(Ji = .2 (10 per cent distance errors) when computing con¬ 
tours in . As the mean may be dominated by a few 
particles with large distances errors that have settled on the 
wrong side of a triple-value zone, we exclude the ten highest 
per-particle in the computation of the mean. Excluding 
the tail of the distribution was found to affect the magnitude 
but not the location of the minimum x^ ■ 

Since our goal is the comparison of predicted and actual 
halo orbits, it was useful to define two further quantitative 
measures of the quality of the orbits, the first, A0, the “di¬ 
rectional error”, being defined as the per-particle average 
angle between the predicted and actual halo velocity vec¬ 
tors at z = 0, and the second. Ad, the “initial displacement 
error”, being the average per-particle distance between the 
actual and predicted halo positions at the first time step. 


3.1 Halo mass assignment: background mass and 
linear growth factor 

While each halo has a well-defined mass according to the 
FOF algorithm, the total mass in haloes is only 37 per cent of 
the total mass in the simulation catalogue. If each halo were 
simply assigned its bare FOF mass, then the mean density 
in the catalogue region as input to NAM would be too low, 
and in the NAM reconstruction the region would behave 
dynamically like a local void. Since the boundary condition 
in the action is the redshift, and since particles at a given 
distance would have higher redshifts in the reconstruction 
owing to the repulsive effect of the void, NAM would predict 
halo distances that are too small. 

In observational catalogues, where the galaxy luminosi¬ 
ties are given and where Dm is assumed to be unknown, 
the standard approach is to multiply each galaxy by some 
constant factor, its mass to light ratio M/L, which may be a 
function of galaxy type. The minimum in x^ as a function of 
M/L will then be governed by a combination of two factors: 
the degree to which the total mass in the catalogue region 
matches the background density set by Dm, and the gravita¬ 
tional dynamics of galaxies, groups and clusters. In practice 
the dominant contribution to y^ is the former factor, which 
is largely independent of the dynamics of the haloes them¬ 
selves. 

In an A-body simulation we define the “halo mass mul- 
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tiplier” Al, the analogue of M/L in the observational cata¬ 
logues, as the global factor by which we multiply the halo 
masses, which in our simulation gives a mass density of 
flm = -3 when M = 2.7. Since Om is a known quantity, 
and since it is also known whether the region under study 
is close to the global mean density, it is possible in this case 
to employ NAM as a probe not only of the total mass but 
also of the relative distribution of mass in haloes relative to 
the background, i.e. the bias. We do this by simply adding 
a smooth constant density component to the value of 
which is input to NAM so that, for any choice of M, the 
total mass density in the catalogue region is matched to the 
input flm. When the total dark matter density in the cat¬ 
alogue region is close to the global mass density, as it is in 
our simulation, this is equivalent to multiplying flm as it ap¬ 
pears in eq. A15 by the factor A4/2.7. When M = 2.7 the 
smooth component vanishes and all of the mass is located 
in the haloes as before. With this reassignment of the cata¬ 
logue mass, a search for a minimum in in the space of M 
becomes a test of the linear bias relation between the un¬ 
derlying mass density and the halo distribution. To reiterate 
this important point: By assuming prior knowledge of flm, 
we can use NAM to probe the relative distribution of dark 
matter in the haloes vs. the background. We will motivate 
this further in the following subsection. 

A further modification to the halo masses is demanded 
by the absence of an analytical zero velocity constraint at 
the initial time step. Since in addition the haloes are mod¬ 
eled as non-evolving point masses throughout their histories, 
NAM solutions involving close interactions between haloes 
at high redshift are commonly seen. While we have not found 
a way, within the framework of the matrix inversion imple¬ 
mentation of the action, to write down an analytical initial¬ 
time zero-velocity constraint, we have found that initial ve¬ 
locities are effectively suppressed by the procedure, already 
employed in modeling the tidal field, of rescaling the halo 
masses in time according to linear perturbation theory. We 
thus multiply the halo masses at each time step by the stan¬ 
dard linear growth factor (see Peebles 1980, eq. 11.15), which 
is zero at t = 0 and unity at 2 = 0. The value of in eq. 
A15 is further rescaled by the same factor to maintain parity 
at each time step with the halo mass density. This ad hoc 
procedure forces the initial velocities to vanish, eliminates 
hard interactions at early times, and significantly reduces 
the nonuniqueness in the solutions, but at the cost of pro¬ 
ducing halo orbits with typical total displacements some 20 
per cent shorter than the actual paths in the simulation. We 
conclude that there is room for improvement in our method 
of defining halo masses, and anticipate that some variation 
on the above theme, or the discovery of an analytical zero- 
velocity constraint, will open the way to further improve¬ 
ments in the reconstructions presented below. 


3.2 Biasing and the mass to light ratio 

Consider a discrete distribution of objects (haloes) sampling 
the underlying mass density field in a region of volume Vb in 
the universe. Let n(x) be the number density of haloes and 
pm(x) be the mass density field, at comoving position x. De¬ 
note the mean values of n(x) and pm(x) inside Vb by n and 


Pm, respectively. Further, choose Vb to be large enough such 
that Pm i® very close to the universal value flmPc. Assume 
that the halo and DM density fields are related by means 
of a linear bias relation between the the density contrasts 
bh(x) = n(x)/n — 1 and (5m(x) = pm(x)/pm “ 1- We write 
this biasing relation as 

(5h = Mm . (2) 

We will express the bias factor b in terms of the assumed 
Ai of haloes as it is invoked in the formalism of NAM pre¬ 
sented here. This is useful in connecting this formalism to 
standard applications of the linear theory of gravitational 
instability as well as to other implementations of NAM (e.g. 
Nusser & Branchini 2000). A natural way to proceed is to 
use the Poisson equation and the bias relation (2) to write 
the peculiar acceleration, g, of an object in terms of bh and 
to compare that with the corresponding expression derived 
from NAM. The Poisson equation is 

V • g = -dTrCpm^m , (3) 

which upon integration and substituting bh = bSm gives 

g(x) = -Gpj3~^ J d®a;'bh(x')K(x, x') (4) 

where K(x,x') = (x — x')/|x — x'|®. In discrete form, 
g(x) = -G^h~^ ^ K(x, Xi) -b ^GpmM^x (5) 

i 

NAM gives, when a = 1, 

ga = -G ^ miK(x, Xi) + ^GPmcX (6) 

i 

where rm is the mass assigned to object i and Pmc = 
V ; rrii/Vo is the mean mass density estimated from the ob- 
jets. Define moi such that for rm = moi we have Pmc = Pm- 
In this case, a comparison between (5) and (6) implies 
moi = Pjnl^ S'Hd 5=1. Consider now Pmc 7^ Pm- The differ¬ 
ence between the mean mass densities must be attributed to 
a uniformly distributed mass component which can be asso¬ 
ciated with low (luminosity) objects that are not included 
in the catalogue. Low mass objects are typically less biased 
than higher mass ones (Mo & White 1996). Comparison be¬ 
tween the uniform terms in (5) and (6) implies 

b = , and rrii = b~^moi . (7) 

Pmc 

In most of the applications of NAM to real data, the 
masses rrii are estimated from the luminosities, Li, by rrii = 
AiLi where A4 is a global mass to light ratio. Assume that 
Pmc = Pm is obtained for Ai = AIo- Then, for a biased 
distribution of objects, the relations (7) imply that the same 
peculiar acceleration is obtained with the global mass to 
light ratio 

Ai = b-'^Aio ■ (8) 
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In our simulation, since 6 ~ 1 on the scale r >2 /i~^Mpc and 
since Mo = 2.7, we will expect to find the best solutions 
when M — 2.7, i.e. when all the mass is located in the 
haloes. 


4 RESULTS 

We searched for a minimum in in fh® two-dimensional 
parameter space of flm and M, holding Hq constant at 70, 
by finding NAM solutions on a grid of 56 points and in¬ 
terpolating a contour map in ■ The results are shown in 
Fig. 4. At 1(7 the best solutions are found at .27 < flm < .45 
and 1.7 < M < 3.9, consistent with the simulation param¬ 
eter flm = .3 and the predicted M = 2.7 (since 6 ~ 1). A 
second set of solutions was computed without the external 
tidal field, and its absence had little effect on the quality of 
the solutions or the location of the minimum in x^- In nn 
idealized scenario where the exact halo masses are known, 
and where the region being reconstructed is a fair sample of 
the global mean, NAM can thus correctly recover both flm 
and b but with an uncertainty of about 50 per cent in both 
values. In an actual observational situation with a compa¬ 
rable number of particles, where significant uncertainties in 
the galaxy -|- halo masses must additionally be taken into ac¬ 
count, it is expected that the uncertainties in the predicted 
values of flm and b will be even larger. It should be remarked 
that the tightness of these constraints is partially a function 
of the size of the catalogue, since sigma in a standard x^ 
measure is an inverse square-root function of the number of 
particles entering into the computation of x^- 

A series of closer looks at the best obtainable NAM 
solution (Figs. 5 through 9), all plotted from the same so¬ 
lution with input parameters flm = .3, M = 2.7, gives a 
good feel for the strengths and limitations of the method. 
Note that for this solution all of the mass is located within 
the haloes. Fig. 5 shows in x-y projection the final positions 
of the haloes in the high-density region of the simulation 
(open circles, with radii proportional to the mass) and the 
error in the predicted positions (line segments). Note that 
the distance errors are purely radial due to the boundary 
condition in the action. The reconstruction tends to be less 
accurate in the vicinity of massive haloes, which is due to 
the presence of triple-value regions where NAM is easily con¬ 
fused. This can be seen more clearly in Fig. 6, which shows a 
view of the “sky” from the reference halo, highlighting those 
haloes with the largest x^- Some of these haloes are found 
near the edge of the catalog region, where interactions with 
particles in the low-resolution region, which is only approx¬ 
imately modelled by NAM, may disrupt the reconstruction. 
The majority of the remaining haloes with high x^ often 
found within a few degrees of the line of sight to a mas¬ 
sive halo; these have typically settled on the wrong side of 
a triple-value region. Haloes with high x^ also tend to be 
found close to each other, suggesting that if the dynamics of 
one influential halo are not correctly modelled it may disrupt 
the accurate reconstruction of several other nearby haloes. 
Finally, a few small haloes with large distances errors can 
be seen in relative isolation, indicating that the interaction 
of haloes with unassociated dark matter particles, and the 



Density parameter 

Figure 4. Mean x^ contours showing the sensitivity of NAM to 

and, through M, to the linear bias b. Solid contours lines 
are from solutions including the external tidal field and dotted 
contour lines are from a parallel set of solutions computed with¬ 
out the tidal field. Contour levels in both cases are .06 or la. 
The lowest computed mean of 1.3 (including the tidal field) 
was found at Om = .3, M = 2.7 (corresponding to b = 1). At 
each grid point a set of 16 solutions was computed, half with and 
half without the tidal field. For values of M different from 2.7 a 
smooth component is added to ensure that the total mass density 
in the catalog matches the background density. The filled square 
at center shows the global simulation parameters, assuming a lin¬ 
ear bias of unity, and the nearby filled circle the parameters in 
the catalog region, which is slightly overdense due to mass inflow 
across the catalog boundary. 

complex formation and merger history of the haloes them¬ 
selves, may prevent accurate NAM reconstruction even in 
the absence of more obvious disruptive factors. 

The top panel of Fig. 7 shows the error in the predicted 
distances as a function of the distance from the reference 
halo. Here the excellent overall prediction of halo distances, 
with typical errors of less than 3 per cent, can be seen most 
clearly. Any mismatch in the halo mass density relative to 
Hm would be revealed here by a tilt in the distance errors as 
a function of distance from the reference halo (with a posi¬ 
tive slope indicating an overdensity and a negative slope an 
underdensity). As the distance errors trace a line with van¬ 
ishing slope, this confirms that the total mass for a choice 
of AI = 2.7 is well matched to Hm. The bottom panel of 
Fig. 7 compares the distance errors with those obtained as¬ 
suming zero peculiar velocities and Hubble-flow distances 
(d = cZiHo). In the latter case the average distance errors 
are 5 per cent. The difference in the sharpness of the peaks 
in the two histograms gives an indication of the ability of 
NAM to correctly model the interparticle dynamics. 

Fig. 8 compares the NAM-reconstructed halo orbits 
with the actual halo orbits, the latter being defined pre¬ 
viously as the CM motion of the dark matter particles com- 
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Figure 5. catalogue vs. predicted halo positions at 2 : = 0 in x- 
y projection, for the best case solution at Qm = .3, A4 = 3.2. 
The dotted outer circle (in green) marks the nominal 26 h~^Mpc 
boundary of the high-density region of the simulation. Actual halo 
positions at the final time step are marked by circles with radii 
proportional to . The free end of the line segments mark the 
positions predicted by NAM (red when the predicted distances 
are larger and blue when they are smaller); the length of the line 
segment thus reveals the radial distance error. 


prising the halo at 2 = 0. The reconstruction is more accu¬ 
rate for the heaviest haloes (top panel) than for the rest (bot¬ 
tom panel): for the former the directional error Ad = 41°, 
while for the latter Ad — 48°. By comparison, the chance 
that a random orbit will have a directional error of 48° or 
less, that is, that a point placed at random within a sphere 
will fall within the volume of the cone swept out by an open¬ 
ing angle of 2 >i= 48°, is about 17 per cent. According to the 
directional error the quality of the orbit reconstructions, like 
X^, is not a particularly sensitive function of A4: far from 
the minimum at the same flm and A4 = 1, Adi — 50° 
for the entire catalogue. Similarly, the average initial dis¬ 
placement error Adi, was 2.5 Mpc/h for the best solution at 
A4 = 2.7 , while far from the minimum at M = 1 it was 2.9 
Mpc/h. A further indication of the approximate character of 
the predicted orbits is that the magnitude of the initial dis¬ 
placement error is comparable to the total distance travelled 
by the typical halo orbit: The centre of mass of the average 
halo in the simulation travelled 3.2 Mpc/h, while the re¬ 
constructed haloes travelled an average of 2.6 Mpc/h, the 
shorter path lengths in the reconstruction being a feature 
of our halo mass assignment scheme as discussed in section 
3.1. While this error may seem large, the chance that a ran¬ 
dom orbit with a total displacement of 2.6 Mpc/h will end 
up within 2.5 Mpc/h of the actual initial position is only 
16 per cent (this is the volume overlap of two spheres of 
equal radii whose centers are separated by a distance equal 
to their radii). A trial NAM reconstruction without the lin¬ 
ear growth factor, assigning the full halo mass at early times, 



o ^ O 


0 100 200 300 

L (degrees) 

Figure 6. A view of the “sky” from the perspective of the refer¬ 
ence halo near the center of the simulation, highlighting the loca¬ 
tion of haloes with poorly recovered distances in the best NAM 
solution. Point size is proportional to the halo mass. Filled points 
indicate haloes with > 1 (radial distance error greater than 10 
per cent). For haloes located towards the edge of the catalogue 
(> 2AMpc/h) these are indicated by filled squares, and towards 
the center of the catalogue (< 2AMpc/h), by filled circles. 


gave similar values for the late-time measures Adi, 

while the early-time measure Adi ~ 4.9Mpc/h, or about 
twice the error. The average total distance travelled by the 
haloes in these solutions was 6.3 Mpc/h, nearly twice as long 
as the actual halo paths and illustrative of the instability of 
the solutions when the haloes are allowed to keep their full 
masses in the initial time steps. 

Finally, Fig. 9 compares three measures of errors in 
NAM reconstructions: Adi and Adi are plotted on the x and 
y axes, respectively, while the point size is proportional to 
X^- Haloes with > 1 shown as filled circles, x^ is fairly 
well correlated to Adi, indicating that reconstructed orbits 
which begin at positions well removed from their actual 
starting positions in the simulation are likely to end with 
inaccurately modeled distances. Significantly, x^ is poorly 
correlated to Adi. This may have been expected since, in 
the absence of nearby haloes constraining the dynamics, the 
motion of a given halo in the plane of the sky relative to 
the reference halo should be fully degenerate. The extent to 
which this degeneracy is broken and the transverse orbital 
motions at the present epoch are correctly recovered is a 
measure of the sensitivity of NAM to the dynamics between 
haloes. The weak correlation between x^ ^nd Adi is the 
clearest evidence that the predicted halo distances alone are 
not a sufficient discriminator of the quality of reconstructed 
orbits. 
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Distance error in units of 


Figure 7. Top panel : Scatterplot of the error in radial distance 
predictions for the best NAM solution, showing good overall re¬ 
covery of the distances. Point size is proportional to halo mass. 
Bottom panel: Histogram in x = — /i9“* )/<Ji of the same so¬ 

lution. For comparison the dotted line shows the errors obtained 
when Hubble-fiow distances are used in place of 

5 DISCUSSION AND CONCLUSIONS 

We have shown, using a catalogue of over 500 dark mat¬ 
ter haloes derived from a large A-body simulation at the 
scale of the local supercluster, that it is possible with the 
numerical action method to reconstruct the full dynamical 
histories of dark matter haloes given the masses and the 
redshift space coordinates at the present epoch. The recon¬ 
struction is most successful in recovering the halo distances 
at the present epoch, with typical errors of less than 3 per 



x| Mpc/h] 

Figure 8. x-y projection of actual (solid lines) and reconstructed 
orbits (dotted lines) for the heaviest haloes (top panel) and for 
all other haloes (bottom panel). Dotted lines indicate the actual 
halo paths to 2 = 20, while solid lines mark the paths predicted 
by NAM. Actual halo positions are indicated by circles of radii 
proportional to the halo mass. Straight radial dotted line seg¬ 
ments connect these positions to those predicted by NAM, as in 
Fig. 5. Heavier haloes tend to have more accurately reconstructed 
orbits: For the heaviest haloes = 41° and for the remaining 
haloes = 48°. 


cent. Individual orbits paths, including the initial positions 
as well as the direction of motion of the haloes at the present 
epoch, are predicted with less accuracy. By varying the rel¬ 
ative contributions to the total mass from the haloes and 
the background, we have also found a way to use NAM to 
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Figure 9. Scatterplot showing errors in the direction of the re¬ 
constructed final-time velocities, as well as errors in the early and 
late time predicted positions, for each of the 533 particles in the 
best NAM solution. The x-axis shows Adi, the error in the initial 
placement of the haloes at the first timestep. The y-axis shows 
A9i, the error in the direction of the velocity vector at the final 
timestep. Halo orbits which most accurately predict the direction 
of the halo orbit at the present epoch are thus found at the top 
of the graph. Point size is proportional x^. 


directly measure the linear bias when the total mass density 
is known, although with an uncertainty of about 50 per cent. 

Given the dynamical complexity of millions of interact¬ 
ing particles, and the sweeping nature of NAM’s central sim¬ 
plifying assumption that galaxy haloes can be approximated 
as discrete, non-merging point masses thronghout their evo¬ 
lution, it is remarkable how successfully the dynamics of a 
many-body system can be reconstructed on the basis of an 
incomplete catalogue of facts. The successes of NAM as it 
has been implemented here are of course partially offset by 
their weaknesses. Among these is the relatively poor qual¬ 
ity of the reconstruction in the vicinity of massive haloes, 
clearly seen in Fig. 5, which shows a breakdown in the non¬ 
linear regime where NAM, which is itself a fully non-linear 
method, might have potentially offered the most insight. In 
these regions is a good indicator of poorly reconstructed 
orbits, but preliminary attempts to use this information to 
nudge haloes into the correct orbits while not imposing any 
additional formal constraints have so far been unsuccessful. 
A second concern is the inability of NAM in many cases to 
isolate, on the basis of alone, predicted halo orbits which 
are moving in the wrong direction at the present epoch. 
Fig. 9 shows, for example, 34 haloes moving more than 90° 
in the wrong direction but with good distances at the present 
epoch and thus low 

The above innacuracies may in part arise from the de¬ 
tails of our implementation, such as our ad hoc procedure of 
scaling of halo masses according to linear theory, and there 


is certainly room here for improvement. The scale of the cat¬ 
alog is also a factor to consider, and in particlar the density 
of mass tracers within it. This analysis should be repeated 
at the scale of the Local Group, where a larger number of 
mass tracers acting within a smaller volume may better con¬ 
strain the dynamics and permit more accurate orbit recon¬ 
structions. It is also possible that, in dynamical systems of 
this complexity, the angular positions, redshifts and masses 
are by themselves insufficient to lift the degeneracies in the 
halo orbits, and that the full three-dimensional velocities 
at the present epoch will be needed to accomplish this. This 
again is work to be undertaken at Local Group scales, where 
next-generation observations from SIM and GAIA hold out 
the promise of multiple galaxy proper motion measurements 
with which to test the NAM predictions. One related con¬ 
cern is that part of the proper motion data may be needed to 
recover accurate orbits, leaving fewer remaining free param¬ 
eters to assist with the more weighty problem of constraining 
individual galaxy halo masses, although it is possible that 
only one of the two components of the tangential velocity 
will be sufficient to break the orbtial degeneracy. Finally, 
inaccuracies in the NAM predictions are doubtless due at 
least in part to intrinsic limitations of the method and its 
assumptions, although we do not wish to suggest at this 
stage, given the work that remains to be done, that that an 
upper limit on NAM accuracy in orbit reconstruction has 
yet been reached. 

We anticipate that work on NAM in the near term will 
lie principally in two directions. The first is an extension of 
the above analysis, with further improvements in the imple¬ 
mentation, to a high-resolution simulation at the scale of the 
Local Group, where present-day three-dimensional veloci¬ 
ties can provide significant additional dynamical constraints. 
The second is a direct comparison of NAM with other recon¬ 
struction methods, both in real space (e.g, Nusser & Dekel 
1992, Gramman 1993, Groft & Gaztanaga 1998, Frisch et 
al. 2002, Mohayaee et al. 2005) and redshift space (e.g., 
Narayanan & Weinberg 1998, Monaco & Efstathiou 1999, 
Mohayaee & Tully 2005), that help to bridge the present-day 
observations of large-scale structure with the initial condi¬ 
tions prevailing in the early universe. 

We acknowledge the support of the Asher Space Re¬ 
search Institute. We would like to thank Felix Stoehr for 
providing us with the snapshots from his simulation. 
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APPENDIX A: THE COSMOLOGICAL ACTION 


The standard cosmological action in an expanding coordi¬ 
nate system, following Peebles (1989), is 


5 = 


/ —A 

“'mE 


mia{t)^ /dxi 


m 


+-E 


rrurrij 

Xi -Xj 




(Al) 


where Xi(t) are the co-moving coordinates, a(t) is the cos¬ 
mological scale factor, and p{t) is the mean background mass 
density. 


The variation in the action, after integration by parts, 


rto 


5S = 


dt 5:>ii 


d 2 dxi 

— — rma -r- J-miflgi 


dt 


dt 


2r dxi 
rriia ox^- 

dt 


*0 


= 0 , 


(A2) 


where 






+ 


4 

3 


TrGpaxi. 


(A3) 


The usefulness of the action in a cosmological context hinges 
upon the observation that the boundary term in (A2) van¬ 
ishes either when (5xi = 0 or a^dxi/dt = 0. The former 
condition obtains when all three position coordinates at to 
are fixed. The latter condition is met automatically because 
a = 0 at t = 0. 

Since cosmological redshifts can be more accurately 
measured than distances, however, it is desirable to find a 


way to fix the observed redshifts as a boundary condition in 
the action and leave the distances as free parameters, while 
leaving the angular position coordinates unchanged. The 
procedure can be transparently carried out in the Hamil¬ 
tonian formulation of the action, where the independent 
status of generalized coordinates q and momenta p can be 
exploited: 



As described fully in Phelps (2002), the coordinate 
transformation which exchanges the roles of the radial ve¬ 
locities and radial distances is carried out by adding a gen¬ 
erating function to the Hamiltonian of the form 

F = gV- (A5) 


where the superscript r refers to the radial component. 
The modified action is 




Pi ■ dqi - ^dt - Vidt - plql. (A6) 

2m 


All boundary terms are evaluated at t = to, and 


y _ Gmimj 

^ <lij ^ 


6 


(A7) 


The subscripts i and j will refer to the particle index, while 
superscripts will refer to coordinate indices. 

In a homogeneous expanding coordinate system with 
scale factor a{t), the physical position is 


ax 


(A8) 


and the physical momentum is 

max — —h max, (A9) 


where x is the coordinate position and where we have re¬ 
defined the variable p = ma^x. Omitting for the moment 
individual particle subscripts, the action is then 
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(—h max) ■ {adx J- axdt) 


—dt ( 


a ^ ^ m .2 2 , xr 


(P"' , ■ r\ 

-h maox aox 

\ao ) 


p ■ dx + madx ■ dx + l^d^x^dt 
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2ma^ 
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— p X — m{x j aoao. 
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(All) 


After integration by parts, and substituting the stan- 
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dard pressureless form of the acceleration of the cosmologi¬ 
cal expansion, 



-f 


A 

3’ 


(A12) 


the action can be written as 


S = 


p ■ dx ■ 


-a 
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+ V ]dt 


-p X 


2ma^ 
r m(x^)^aodo 


(A13) 


where V is the gravitational potential modihed by a 
smoothly distributed homogeneous background density, 


Vi 
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a rriiXi. 


(A14) 

(A15) 


The cosmological constant is no longer explicitly found 
in the action, as expected since it cannot produce peculiar 
accelerations in co-moving coordinates. Its influence is felt 
in the timescale of the integration. 

The action has thus been put in the standard Hamilto¬ 
nian form, with conjugate variables x and p, and with the 
addition of extra terms outside the integral, the first repre¬ 
senting the generating function of the canonical transforma¬ 
tion, and the second arising from the change of boundary 
conditions (for further details see Phelps 2002). 

To recover the correct equations of motion, these bound¬ 
ary terms must vanish in the variational derivatives, and so 
the constraint which was implicit in the transformed coor¬ 
dinate system - that the final-time velocity of each particle 
relative to the origin is equal to its observed redshift - must 
be imposed by hand in the original coordinates. As the origin 
of the coordinate system is typically taken to be the Milky 
Way, its motion along the line of sight must be subtracted 
out: 


Vpec + Hor - Vo ■ X = z, 
dox^ 

H- Vo ■ X = z, 


(A16) 


lltLtQ QiQ 

where uo is the coordinate velocity of the Milky Way. 

With the equation of constraint (A15), the action is 
written 

o — —aoao[x ) — maozx + maoX vo ■ x 


+ / (p ■ dx — H\\J). 


(A17) 


By construction, solutions found by minimizing this 
modihed action will contain particle orbits whose hnal time- 
step angular positions and radial velocities relative to the 
reference particle will equal the input angular positions and 
redshifts. Distances and transverse velocities at z = 0, as 
well as the orbit history of each particle, emerge as predic¬ 
tions. 

The rest of the problem is in principle a computational 


one, and consists in devising a numerical approach to hnd- 
ing solutions to the variational integral. The hrst numer¬ 
ical applications of the action method used a set of htting 
functions to characterize the orbits, and the action was mini¬ 
mized with respect to the coefficients of these functions. The 
procedure employed here follows Peebles (1995), where the 
positions are independent at each time step, allowing for the 
creation of sharper orbits. 

In this implementation the position coordinate of 
the particle at the time step is written as a;“„. If 
there are N time steps, then a:“jv+i is defined to represent 
the present-time positions (the hnal time step is treated dif¬ 
ferently from the others, due to the boundary conditions). 
Time increments are dehned every half step, so that a = 1 at 
t 2 N+i (“fo”) and a = 0 at the hrst time step. The momen¬ 
tum is then dehned, as in a standard leapfrog integration, 
at the half-steps, as p“ 2 n- With this notation the action in¬ 
tegral becomes a sum and the per-particle action is written 
as 




VHi . / r \2 r 

2 d2N-\-l{^Xi TfliXij\f^iZ 

+miXi jv+lUO iV-l-1 • Xi N+1 
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+ V„{t2n + 1 — t2„-l) 


(A18) 


where the gravitational potential V is dehned as in equation 
(A15). 

The condition that the action be stationary with respect 
to the p-derivative yields an expression for p“ 2 n in terms of 
X and t: 


dSi 

dp?2, 


n+1 


^2 i^Zn+l — t2n-l) 


0, (A19) 
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(A20) 


This is substituted back into the action, which is then 
expressed solely as a function of x: 

Oi — -^o,2N+i(Xi — rriiXi ff^iZ 

+miXi jv+lUO iV-l-1 • Xi N+1 

VlliOj2n {Xi n+1 Xi n) 

2 t2n + l — t2n-l 

— V„(t2n+1 — t2n-l)\ ■ (^21) 

To hnd stationary points in the action through matrix 
inversion, by means of the iterative approach described in 
Peebles (1995), the hrst and second derivatives of the action 
with respect to the position variables at each time step must 
be computed. Calculation of the derivatives is considerably 
simplihed by hxing the orbits of all particles except the one 
whose orbit is being solved, as the off-diagonal terms can be 
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ignored. The derivatives at the Hnal time step are computed 
separately. 

The gradient of the action for all time steps except the 
last is: 

OSi _ 2 is^in ^in—l) 

9 a;“„ _ i 2„_3 

2 {Xi„+1 - xfP) 


d^Si 


mia2pfXi iv+i 
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for the last two time steps. 


(A22) 


and the second derivatives are 
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for different time steps. 

The gradient and second derivatives of the potential are 
BVn Gmitrij x'ji „ 
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At the final time step, when i 1, only the components 
of the gradient in the radial direction are needed: 

dS 

j— -= mi{HoXi — z + Vo iv+r • a;“jv+i) 
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where Xi n+i is referred to the origin and so is purely radial. 
The second derivatives are 
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